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We study the local density of states at and around a substituting impurity and use these results to compute 
current versus bias characteristic curves of Scanning Tunneling Microscopy (STM) experiments done on the 
surface of graphene. This allow us to detect the presence of substituting impurities on graphene. The case of 
vacancies is also analyzed. We find that the shape and magnitude of the STM characteristic curves depend on 
the position of the tip and on the nature of the defect, with the strength of the hinging between the impurity and 
the carbon atoms playing an important role. Also the nature of the last atom of the tip has an influence on the 
shape of the characteristic curve. 

PACS numbers: 73.20.Hb, 73.23.-b, 81.05.Uw 



I. INTRODUCTION 

Graphene''^ consists of a monolayer of carbon atoms form- 
ing a two-dimensional honeycomb lattice. It has been inten- 
sively studied due to its fascinating physical properties'* and 
potential applications. The honeycomb lattice consists of two 
triangular sub-lattices and this is responsible for the linear dis- 
persion of the low-energy excitations and for a pseudospin de- 
gree of freedom for electrons in graphene. Many of the novel 
properties of graphene follow from these two facts. Because 
of the Dirac spectrum, disorder can have a significant effect 
on the electronic properties of graphene, the effect being es- 
pecially strong when the chemical potential crosses the Dirac 
point. Extrinsic disorder in graphene can be in the form of 
impuritieS )'*i^'^'^-^ topological defects^iiMiii^ edgesr^iii and 
substrate corrugations.'^ In addition, there is also disorder in 
the form of intrinsic ripples in the structure of graphene ^'^-'^''^ 
Disorder in graphene occur naturally, but can also be induced 
if this is advantageous, to tailor its transport properties. This 
is the case for the recently produced material graphane.'^ 
Among the several possibilities, the replacement of a carbon 
atom by a different atom can occur. Atomic substitution in a 
carbon honeycomb lattice is chemically possible for boron (B) 
and nitrogen (N) atoms. There have been several experimen- 
tal studies of B and N substitution in highly-oriented pyroly tic 
graphite, graphitic structures^^^ and nanoribbons.^^* 

In a previous publication^! the problem of chemical substi- 
tution in graphene has been considered, and the local density 
of states (LDOS) and local electronic structure and charge dis- 
tribution have been numerically calculated. In this work, we 
extend the calculation of the spatial dependence of the LDOS 
at and around the impurity using analytical methods and ex- 
tending the calculations for energies way beyond the Dirac 
point, an essential ingredient for the calculation of Scanning 
Tunneling Microscopy (STM) currents at finite bias. Using 
these results we perform theoretical calculations of the tunnel- 
ing current versus bias characteristic curves of STM studies 
on graphene surfaces with dilute substitution impurities. The 
STM technique is one of the most powerful tools for study- 
ing surfaces. This is particularly convenient for the study of 
graphene, which is two-dimensional and exposed to the STM 
tip. There has been several STM studies of both graphene 



grown epitaxially on SiC ) ' ' and mechanically exfoliated 
graphene on Si02 i'^'^^'^^'^^'^' Our results can be used to in- 
terprete the STM signal when the tip comes close to a sub- 
stituting impurity. A related study, computing the LDOS for 
point deffects in graphene, using first principle methods, was 
recently performed. A study of the LDOS starting from the 
Dirac equation considering the effect of magnetic impurities 
was also recently considered. The effect of local potentials 
with finite strength was studied in Ref. |34|. 

In Section |ll] we present the tight binding model and the 
Green's function formalism used. For a single localized sub- 
stitutional impurity, the Green's function can be obtained ex- 
actly in closed analytical form. In Section Hill we discuss re- 
sults for the LDOS and in Section |IV] the results for the tun- 
neling current are presented. Section |V] contain further dis- 
cussions and conclusions. 



II. CALCULATION OF THE GREEN'S FUNCTION FOR 
DISORDERED GRAPHENE 

A. Basic definitions 

Let us consider that a carbon atom on the otherwise perfect 
lattice of graphene has been replaced by a different type of 
atom. What will happen is a renormalization of both the on- 
site energy and the hoping parameter between that atom and 
the nearest neighbor carbon atoms. If the hoping parameter of 
clean graphene is t then the hoping between the impurity atom 
and the carbon atoms can be parameterized by an additional 
constant denoted to, as shown in Fig. ([T]i. The value of to 
can be varied: if the impurity atom is strongly coupled to the 
carbon atom, to will be negative; if, on the other hand, the 
coupling is weak to will be positive. The Hamiltonian of the 
system will be that of clean graphene and a term representing 
the renormalization of the hoping, as explained below. 

Let us first introduce some definitions for latter use. The 
honeycomb lattice has a unit cell represented in Fig. [T]by the 
vectors ai and 02, such that \ai \ = |a2| = a, with a ~ 2.461 
A. In this basis any lattice vector R is represented as 

R = nai + ma2 , (1) 
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FIG. 1: (color on line) The honeycomb lattice with a substituting 
atom replacing a carbon atom. The local hopping parameter changes 
from t to t — to around the A atom. We assume all over t = 3 eV. 



with n, m integers. In cartesian coordinates one has, 

ai - Y(3,V3,0), a2 - y(3,-\/3,0), (2) 

where ao — a/\/3 is the carbon-carbon distance. 

If periodic boundary conditions are used, the Bloch states 
are characterized by momentum vectors of the form 



mi 7712 

k — — Oi H 02 , 

^1 A^2 



(3) 



with 771 1 and 7^2 a set of integers running from to iVi — 1 
and from to A^2 ~ 1, respectively. The numbers A^i and N2 
are the number of unit cells along the Oi and 02 directions, 
respectively. The total number of unit cells is, therefore, Nc — 
NiN2- The reciprocal lattice vectors are given by: 

bi = |^(1,V3,0), b2 = |^(1,-V3,0), (4) 

The vectors connecting any A atom to its nearest neighbors 
read: 

= y(-l,\/3,0) = i(ai-2a2), (5) 
^2 = y(-l,-\/3,0) = i(a2-2ai), (6) 



^3 = ao(l,0,0) = -(ai +02) 



(7) 



Using the above definitions the Hamiltonian for this prob- 
lem can be written as 



Hq = -tJ2HR)bHR)+a{R)b\R-a2) 



R 



+ a{R)b\R-ai) + h.c] , 



which represents clean graphene (the spin index is omitted for 
simplicity of writing), and the perturbation due to the impurity 
atom is 



a(0)6t(-ai) 



a(0)6t(- 
~h.c.]. 



-a2) 



(9) 



It is assumed that the impurity atom is in the unit cell R = 0, 
but there is no loss of generalization due to this choice. In 
the particular case of to — t, the scattering term Vi represents 
a vacancy, since the impurity atom is completely decoupled 
from the carbon atoms. It is important to keep in mind that in 
a given unit cell both A and B type of atoms are described by 
the same vector R. 

The calculation of the electronic properties of graphene re- 
quires the calculation of the corresponding Green's functions, 
whose definitions are 

Gaaik,q,T) = -(T[afc(T)at(o)]) , (10) 

Gtbik,q,T) = -(T[6fc(T)fot(o)]) , (11) 

Gabik,q,T) = -(T[afc(T)fot(o)]) , (12) 

Gba{k,q,r) = - (T [6fc(r) ^ (0)] ) , (13) 

and their Fourier transform to the Matsubara representation 
are given in the Appendix lAl for the case of a perfect lattice. 
Of particular interest to our calculations is the momentum in- 
tegral of the retarded diagonal Green's function 



(14) 



The integral is best performed using the density of states of 
the honeycomb lattice. Since we want to take into account 
the non-linearity of the bands, which allows us to describe the 
properties of the system at large energies and not only close to 
the Dirac point, we have to use for the density of states an ex- 
pression that goes beyond the usually used linear dependence 
of this quantity on energy. In a previous work^I we derived 
an expansion for the density of states (per unit cell, per spin) 
valid for energies up to 3 eV, reading (E ~ hio) 



P{E) 



2E 



2E^ 



lOE^ 



The imaginary part of G^^(cj) reads 



(15) 



(16) 



and the real part has the form 

^G%{lu) = Pi{huj) + P2{hLo) In ■ 



where Pi (x) and P2 (x) are polynomial functions given by 
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(17) 



X 



X 



-Dt 



X 



27i4 V2 

3 5 
+ 



X 



(8) 



Dl 2,f^Dl 27£|2t4 ■ 
The energy Dc is a cut-off energy chosen as — \fiT:t^. 



(18) 
(19) 
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B. Exact Green's Functions This is best accomplished using the equation of motion 

method. The equations of motion for the Green's functions 
We now want to determine the exact expressions for the ^an be readily established, and read 
Green's functions in the presence of the substituting atom. 



iuJnGAA{l-On,k,p) = 5k,p + t(t){k)G BA^I^n, ^iP) ^ ^^^^ (j){k' )G BAji^n, k' , p) , 

iu}nGBA{t^n,k,p) = t(j}* {k)G AA{^n-,k,p) 

iuj.,iGAB{i^n,k,p) = t(f){k)GBB{i^n,k,p) - -^^(l){k')GBB{(^n,k',p) , 

iuJnGBB{(^n,k,p) = 



to 



5k,p + t^i* {k)GAB fc, p) 
I 



Nr. 



{k)^GAB{^n,k' ,p) . 



(20) 
(21) 
(22) 
(23) 



The complex number 4>{k) is defined as 



m 



l + e 



ik-ai _j_ ik-a2 



(24) 



which is the form factor of the three B atoms as seen by an 
atom in A. The above set of equations can be solved exactly. 
The fact that the scattering term Vi depends on 4>{k) and that 
a single impurity breaks particle-hole symmetry of the system 
implies a complex form for T-matrix. In fact, the general ex- 
pression for the Green's functions does not have exactly the 
same form as in the case of one-band electrons. The solution 
of the equations of motion is rather lengthy and in the course 
of the solution we use the two following identities: 



k,p) = l + {t-to)Y,<t>ik)GBA {uJn, k,p) . 
k 

(25) 



k 

and 



toz = -[{iujnf + zto{t-to)]^(l}ik)GBA{^n,k,p) 

k 

+ iu;„t^|0(fe)pGAA(w„,fe,p), (26) 

k 

with 

^ k 

The use of the relations (|25l l and 

^GAA{^n,k,p 



leads to 

Ni{p,uj„] 

D{LOn) 



with the following definitions 

Ni{p,UJn) = {t-to)G%{uJn,p)+toG%{uJn). 



(28) 



(29) 



i^K) = {t- io)(l - to/t) + iu;„i2to - tl/t)G%{LUn)(?0) 

G"AAM^^J2^AAi^",k), 



(31) 



and 



G%K) - 1^ E \m\'G%iu;n,k) . (32) 

^ k 

The result ( |28] ) has the appropriate limiting behavior: when 
to t one obtains {iujn)~^, which agrees with ( |25] ): when 
to ^ one obtains G'^j^{ujn,p), which is the result for the 
perfect lattice. 

Finally, the exact solution for Gaa{^ji, k,p), considering 
an arbitrary value of to, has the following structure 

GAAi^^n, k,p) = Sk^pGAAi^n, k) + G% (cj„ , fc)S(cj„) 
+ g{uJn) + GAA{uJn,k)TAA{k,P,UJn)GAAi^n,P) , (33) 

with 

1 toGAAi^n) 



E(w„) 



Nc t D{0Jn) 
1 toil -to/t) 



and TAA{k, p, cj„) given by 

^ „ ^ 1 tjto - tmk)\' - i^^^n)' 
TAA{k,p,UJn) = —to : pr^ — T 



(34) 



(35) 



(36) 



Clearly, the solution ( [33l l does not have the simple form of 
the T— matrix as in the case of the solution of the scattering 
problem by a local potential in the single energy-band case. 
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The calculation for GBsi'-^n, k,p) proceeds along the same 
lines and gives 

+ G^BK,fc)ri3s(fc,p,w„)G%K,p), (37) 
with TBB(fc,p, w„) given by 

1 r{k){2t - to)tto(t>{p) 



TBB{k,p,UJn) 



Nr. 



iu!nD{LLln) 



(38) 



It is interesting to note that is the case of Eq. (|37| i the 
T— matrix has the traditional form. So, as long as one does 
not sit on top of the impurity atom, the scattering equations 
are similar to those of the scattering of the electron gas by a 
local potential. The exact results (|33] ) and ([37]) is what is need 
to discuss the behavior of current versus bias in STM exper- 
iments when the microscope tip comes near to a substituting 
defect. 



III. LOCAL DENSITY OF STATES 

As we will see later, the properties of the STM current will 
depend on the local density of states below the tip of the mi- 
croscope. We therefore have to compute this quantity for dif- 
ferent circumstances. Because each unit cell is identified by 
a single vector R, the local density of states (per spin) at the 
atoms A and B of the unit cell localized in the position R is 
defined as 



p,j;{R,Uj) = 



1 



ttN, 



■ImGxxiR, R,^) : 



(39) 



with Gxx {R, R, ^) {x = A, B) obtained from 



Gxx{R.R,yJn) = ^e'('^-P)-^C(fc,p,c^„), (40) 

after the usual analytical continuation tj„ — > a; + «0+ of the 
Matsubara Green's function. 

Let us first consider the case — t- With this choice the A 
atom is disconnected from the rest of the lattice and in this sit- 
uation the momentum sum over the full GAAi^^n, k, p) reads 



1 



(41) 



which corresponds to an isolated atom. As a consequence the 
local density of states is a delta function and no charge trans- 
port can take place through that atom. In the material this 
situation never happens and therefore we can interpreted this 
result as the case where there is a vacancy in the lattice. There- 
fore the presence of a vacancy can be detected by its influence 
on the electronic density of states of the neighboring atoms 
(the B atom in the case of Fig. [Til. Using only the first term in 
Eq. dTsl l it is possible to derive a relatively simple expression 
for the local density of states at the B atom near the vacancy, 
reading 



Pb(0, w) 




huj 



with 



1 + ^ In^ 



hhJ 



t 



(42) 



(43) 



It is also possible to determine general analytical expressions 
for /Oa(0, lo) and /Os(0, oj) given an arbitrary value of t^, but 
due to the form of D{iujn) which in this case depends on both 
t and to, the final result is somewhat cumbersome. However, 
if really needed, it is straightforward to use the full equations 
for Gaa and Gbb and the given expressions for Gaa{^^) to 
write down analytical expression for the density of states. 
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FIG. 2: (color on line) Local density of states (LDOS) at the unit cell 
R = Qiov the A and B atoms and different values of to in electron- 
volt. In panel (a) we plot the density of states of clean graphene and 
the LDOS at the B atom near a vacancy at the A atom of the same 
unit cell. For to = 1, Pb{0, u>) is multiplied by 10 for clarity. 

In Figure [21 we plot the local density of states for different 
values of to. In panel (a) of that figure we represent the density 
of states at a S site in the unit cell where a vacancy exists in 
the A atom of the same unit cell. The most significant feature 
is a development of a logarithmic divergence at zero energy. 
Therefore the clean density of states (also shown for compari- 
son) is strongly modified by the presence of such a strong po- 
tential. For a moderate value of to the change of the hoping is 
not strong and the density of states at both the A and B atoms 
retain the linearity of the clean density of states close to the 
Dirac point. However, the absolute value of the two density of 
states are not the same and a clear deviation from linearity is 
seen to take place for lower energies when compared with the 
clean case. If the impurity atom binds strongly to the carbon 
atoms (to = — 1) there is an increase of the density of states at 
the neighbor atoms at the expenses of the density of states of 
the impurity. On the other hand, if the impurity binds weakly 
to the carbon atoms (to — 1), the density of states increases 
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at the impurity and strong resonances develop on the impu- 
rity density of states. Clearly these different behaviors of the 
density of states will show up in the tunneling current. 

It is also instructive to compute the local density of states 
as one moves away from the defect. This amounts to perform 
the Fourier transforms for finite R in Eq. ( l40l i. In the case 
of Gaa{R, R, ^) the calculation of the integrals is facilitated 
by the fact that the momentum space Green's function on the 
sub-lattice A depends on |(/'(fc)| only, whereas in the case of 
the Green's function on the B sub-lattice this in not the case. 
Besides for finite R the useful relation iA9i is valid no more. 
However, for large R — \R\ values when compared to a, i. e. 
a way from the defect, the details of the lattice are no longer 
important and the behavior of the local density of states at 
both the A and B sub-lattices must be similar We will there- 
fore give the density of states as function of \R\ for the A 
sub-lattice only. Performing the Fourier transforms using the 
Dirac cone approximation (this is not a restriction, but helps 
keeping the final result not too cumbersome) one obtains (for 
the retarded function) 



G,,4R,R,Lo)^NcG%iio) + N, 



Qiio) , (44) 



where vp — 3aot/{2h) and Q{uj) can be written as 
where I{lu) is defined as 



I{uj) = /o(w) - i7rsgn(cj)Jo(|w|i?/t;i^) 



(45) 



(46) 



with Jo (a;) the Bessel function of integer order n = 0, and 
finally /o(ijj) is the Cauchy principal value of the integral 



Jq - 



(47) 



which can be evaluated using standard numerical methods, 

and we also have fee = 2ypK j {\J Z^f^ao), and a = ujR/vf- 
In Figure[3]we depict the local density of states at the carbon 
atoms located in the A sub-lattice, computed using Eq. ( l39l l. 
The typical oscillations in the density of states close to im- 
purity centers are clearly seen. Note that as one moves away 
from the impurity the density of states approaches that of the 
clean system. The Friedel oscillations in graphene can be ob- 
tained from integrating the density of states up to the Fermi 
energy. 



IV. CALCULATION OF THE TUNNELING CURRENT 

We now want to compute the tunneling current between 
the tip of an STM microscope and graphene, when the tip is 
closed to an impurity atom. We model the tip by a one di- 
mensional tight-binding system, a standard approach. ^^'^^ Our 
results will not depend significantly from this choice as long 
as we assume a large value for the hoping parameter of the 
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FIG. 3: (color on line) Local density of states (LDOS) at distance 
R from unit cell J? = 0, where the impurity is located for different 
values of to. The energy is htu = 0.5 eV. The case of the vacancy is 
also illustrated. 



electrons in the tip. This choice essentially corresponds to 
represent the tip by a metal with a large bandwidth. Given the 
above, the Hamiltonian for the tip has the form 

-1 

H=-V [c\n-\)c{n) + c{n)c\n~l)]+Ho, (48) 

n— — oo 

and i/o representing the Hamiltonian of the last atom of the 
tip. This choice corresponds to the assumption that the surface 
atom of the tip has a different nature from the atoms in the bulk 
of tip. We therefore represent Hq by 

Ho = eoc^(0)c(0) - M^i[ct(-l)c(0) + c{Q)c\-l)] . (49) 

It is essential that the Hamiltonian of the tip to be represented 
by a semi-infinite metal, since otherwise there would be elec- 
trons reflected at the far end of the tip modifying in this way 
the tunneling current. Finally we need to include the tunnel- 
ing of the electrons of the tip to graphene. There is a number 
of ways we can do this. Here we assume that the coupling is 
made directly either to the impurity atom or to the next neigh- 
bor carbon atom. This choice corresponds to probing the local 
electronic properties at or around the impurity. More general 
types of coupling are easily included in the formalism. We 
write this coupling as 



-W2[c\0)d{Q) + d\Q)c{Q)] ., 



(50) 



where the operator d{Q) can represent either the impurity atom 
at the A sub-lattice or the carbon atom at the B sub-lattice. 

Since the Hamiltonian of the problem is bilinear we can 
write it in matrix form (of infinite dimension) as 



H = 



Hb Vl 
Vl Ho 
Vb. H„ 



(51) 



where the matrices Vl and Vr represent the coupling of the 
last atom in the tip of the STM microscope to the bulk of the 
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tip and to graphene, respectively, and Hb and Hg stand for the 
bulk Hamiltonians of the tip and of graphene, including the 
impurity potential (|9]l, respectively. 

The tunneling is a local property, controlled by the coupling 
of the last atom of the tip to the bulk atoms and to graphene. 
Since we want to compute local quantities, this is best accom- 
plished using Green's functions in real space. The full Green's 
function of the system is defined by 



{lE + iQ+ -H)G^ 



(52) 



where we have chosen the retarded function (denoted with the 
+ superscript), and 1 is the identity matrix. The matrix form 
of the Green's function is 



Gbb GbO Gbg 

Gob Gqo Gog 

Ggb GgQ Ggg 



(53) 



The quantity of interest is Gqo, 
the form 



which can be shown to have 



G+ =iE + iO+ - eo - - 



(54) 



where the matrices and E]j are the self energies and have 
the form 



S+ = W^Gt 



(55) 



where the Green's functions G+ and G+j, are the surface 
Green's function of the Hamiltonians Hb and Hg at the impu- 
rity unit cell {x = A, B), respectively. Note that the quantity 
G+j. is computed using Eq. ( |44| ) making R — 0. It is possible 
to find a close form for G^^, as we have shown in the previous 
section. Also for Gf a close form exists^2i^2il£ 



G+ = [E + iO+ -V^Gtr^ . 
The solution of Eq. ( |56] l is elementary and reads 



G a 



21/2 2^2^ 



for E^ < 41/2 and 



G ^ 



E 

2y2 



sgn(^) 
2y2 



^£:2 _ 41/2 



(56) 



(57) 



(58) 



for £;2 > 4y2. 

Our goal is to study the STM current at finite bias, which is 
a particular case of non-equilibrium transport. This is done us- 
ing the non-equilibrium Green's function method or Keldysh 
method. This method is particularly suited to study the regime 
where the system has a strong departure from equilibrium, 
such as when the bias potential VJ, is large. We consider, how- 
ever, that the system is in the steady state. Since the seminal 
paper of Caroli et al. on non-equilibrium quantum transport,"^' 
that the method of non-equilibrium Green's functions started 
to be generalized to the calculation of transport quantities of 
nanostructures. There are many places where one can find a 
description of the method,"^^!"*^ but a recent and elegant one 



was introduced in the context of transport through systems 
that have bound states, showing that the problem can be re- 
duced to the solution of kind of quantum Langevin equation.— 
The general idea in this method is that two perfect leads 
are coupled to our system, which is usually called the device. 
In our case the device is defined by the last atom at the tip 
of the microscope. The Green's function of the device has 
to be computed in the presence of the bulk of the tip and of 
graphene. This corresponds to our G^g Green's function. Be- 
sides the Green's function we need the effective coupling be- 
tween the last atom of the tip and the bulk atoms as well as 
that to the graphene atoms, which are determined in terms of 
the self-energies 



^L/R — 



— fs+ 



- S 



L/R 



)• 



(59) 



Therefore the effective coupling depends on the surface 
Green's function of the tip and of graphene. According to the 
general theory, the two systems (bulk of the tip and graphene) 
are in thermal equilibrium at temperatures and chemical 
potential fii^ and are connected to the system at some time 
tfl. The bottom line is that the total current through the device 
is given by (both spins included) 

2e f°° 

J=-y dET{E)[f{E,^lL,TI:)- f{E,^lR,TR)], 

(60) 

where f{x) is the Fermi-Dirac distribution and the transmis- 
sion T{E) is given by 



T(i?)=4^2Tr[r^G+rflGJ. 



(61) 



Performing the trace (which in this case is only a product of 
complex numbers) we obtain 



T{E) 



m^wi^Gt^G+jG+^\ 



(62) 



Figure|4]represents T{E) in several conditions. In the clean 
case, we see the effect the Dirac point has on the transmis- 
sion, leading to a suppression of the tunnelling for i? ~ 0. It 
is also clear that the parameters characterizing the last atom 
of the tip strongly influences the form of T{E), leading to an 
asymmetry between negative and positive energies due to the 
finiteness of eq- When a vacancy is present at the A site, the 
tunneling through the nearest B sites is strongly modified rel- 
atively to the clean case, with strong tunneling taking place 
at energies very close to the Dirac point. When an impurity 
atom sits at the A site the behavior of the tunneling probabil- 
ity depends on the coupling between the impurity and the car- 
bon atoms. In the case of weak coupling to the carbon atoms 
(to > 0) the tunneling through the A atom is faciUtated rela- 
tive to the tunneling between the nearest carbon atoms. When 
the impurity is strongly coupled to the carbon atoms (io < 0), 
the reverse happens. This behavior is easily understood, since 
strong (weak) coupling to the carbon atoms is equivalent to 
a weak (strong) coupling to the tip of the microscope, lead- 
ing to weaker (stronger) tunnelling probability from the tip to 
graphene. 
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FIG. 4: (color on line) Transmission coefficient T{E) as function of 
the energy for zero bias. The fixed parameters are: t — 3.0, V — 1, 
Wi — 0.9, 14^2 = 0.2, all in electron-volt. Panel (a) depicts the case 
of clean graphene and two different values of eo. Panel (b) depicts 
T{E) through a B site when a vacancy sits at the A site. Panel (c) 
depicts T{E) through A and B sites, when an impurity sits at the A 
site, taking to = 1 eV and eo = 0. Panel (d) is the same as panel (c) 
for to = -leV. 

All the above discussion applies to zero bias voltage. When 
finite bias voltage is applied between the tip and graphene, 
the behavior of the curves change. We will consider that 
the chemical potentials of the tip and of graphene differ by 
the electrostatic energy eVb, where Vb is bias potential. This 
amounts to change the on-site energies relative to their value 
in equilibrium. We choose to change the on-site energies of 
the tip by eVb/2 and those of the carbon atoms by —eVb/2. In 
Figure|5]we depict T{E) for a finite value of the bias Vb- The 
most distinctive difference relative to the case of zero bias is 
the energy asymmetry induced by the bias relatively the case 
of zero bias, even when eq = 0. 

The calculation of the current also depends on the chemical 
potential. In the case of graphene this can be tuned by a back 
gate Vg. The relation between the back gate voltage and the 
chemical potential of graphene obeys the relation 

H = VFh\l^^^^ , (63) 

which is numerically equal to ~ O.OSv^g in electron-volt. 
The parameters in Eq. ( |63] ) are e = 3.9, the dielectric constant 
of Si02, d ~ 300 nm, the thickness of the Si02 substrate, and 
e the elementary charge. Taking a typical value of Vg = 100 
V we obtain /i — 0.3 eV, which is the value we assume for the 
chemical potential in the calculations below. Since we will 
perform our calculations at zero temperature, the cuiTent is 
given by 

J = -r / T{E, Vb) . (64) 

" Jfj,+eVi/2 

The results for the calculation of J versus Vb is depicted in 
Fig. |6] Looking at them, we see that there is an asymmetry 



FIG. 5: (color on line) Transmission coefficient T{E) as function 
of the energy for finite bias Vb = 0.5 volt. The parameters and the 
description are the same listed in the caption of Fig. |4] The energy 
scale in the case of the vacancy was increased in order to see the low 
energy behavior close to the Dirac point. 

between the negative and positive values of Vf,. For the clean 
case, it is clear that the magnitude of the J current does not 
depend much on the chemical potential, and therefore on the 
gate voltage Vg. The same is true for the vacancy case. Com- 
paring the curves for the vacancy and for the clean case, we 
see that the order of magnitude of the current is the same, but 
the shape of the curves is notoriously different from that of the 
clean case. When the impurity is weekly coupled (to > 0) to 
the carbon atoms, the tunneling current is facilitated through 
the impurity atom, as can be seen from panel (c) of Fig. |6] 
with an absolute value five times larger for Vf,= l volt. Also a 
clear jump is seen in J (panel (c)), a finger print of the res- 
onances in the density of states seen in panel (c) of Fig. |3] 
One note that in Fig. |3]the resonances are symmetrically po- 
sitioned relatively to the Dirac point, but this is not so for the 
jumps in the current. The reason is due to finite on-site energy 
at the atom in the tip of the microscope. On the other hand 
for strong coupling, panel (d) of Fig. |6] there is not much dif- 
ference from the clean case, apart from the fact that tunneling 
cuiTent through the impurity or the neighbor carbon atom has 
different magnitude. 

Another important quantity to fully characterize the STM 
current is the shot noise 4^ For interacting systems this quan- 
tity gives information on the possible existence of quasi- 
particles with fractional charge. On disordered systems with 
no interactions, information on transport open channels can 
be obtained. For fermionic non-interacting electrons at zero 
temperature the shot noise is defined as"*^ 

The relevant quantity is not S directly but the Fano factor- 
defined as 
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FIG. 6: (color on line) Current J as function of the bias voltage 
V(,. The parameters are the same listed in the caption of Fig. |4]and 
£() = 0.2 eV. Panel (a) represents the clean the case for two values 
of the chemical potential /.i. The case of the vacancy is represented 
in panel (b), also for two values of the chemical potential. Panel (c) 
represents the case to = 1 eV for = 0.3 eV, and panel (d) is the 
same as (c) for % = —1 eV. In panels (c) and (d), A and B refer to 
tunneling through the atom siting on sub-lattice AoxB, respectively. 



When the transmission TiE) is strongly reduced we have 
— > 1, and noise is said poissonian. On the other hand, if 
the system has a finite density of open channels, T{E) 1, 
we have F < 1 due to [1 - T{E)\ < 1. In Figure |7] we plot 
the Fano factor as function of the bias voltage. 
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FIG. 7: (color on line) Fano factor F as function of the bias voltage 
Vf,. The parameters are the same listed in the caption of Fig. |4]and 
eo = 0.2 eV. Panel (a) represents the clean case for two values of the 
chemical potential /i =0 and /i —0.3 eV. The case of the vacancy is 
represented in panel (b), also for two values of the chemical potential. 
Panel (c) represents the case to = 1 eV for /i — 0.3 eV, and panel 
(d) is the same as (c) for to — ~1 eV. In panels (c) and (d), A and 
B refer to tunneling through the atom siting on sub-lattice A or B, 
respectively. 

Clearly, for all case but the vacancy, the Fano factor is close 
to the poissonian values F — 1. When Vb ^ and the chem- 



ical potential is zero, the value of T{E) quite small due to the 
presence of the Dirac point, and the Fano factor take the lim- 
iting value F = 1. For the vacancy, we obtain much smaller 
values of F, which is indication of the presence of the strong 
resonance seen in the local density of states for this case, 
which leads to an increase of the transmission. It should be 
noted that in the cases of finite to, the concavity of the curves 
for F as function of Vb allows to distinguish the case of week 
coupling to > from the case of strong coupling to < 0. 



V. DISCUSSION AND CONCLUSIONS 

We have studied in detail the electronic local density of 
states in graphene close to a substituting atom. This quan- 
tity turns out to be important in STM experiments, since the 
tunneling of the electrons from tip of the microscope to the 
sample is determined by this quantity. The shape and mag- 
nitude of the tunneling current depends on the nature of the 
substituting atom. In the case of the vacancy the shape of 
the tunneling current has a very different signature from the 
other cases. Also the magnitude of the current depends on the 
strength of the bonding between the impurity and the carbon 
atom. 

The curves we have presented here are only indicative of the 
general behavior of the density of states and of the tunneling 
current. In order to obtain experimentally relevant curves we 
should have real numbers for the parameters of the tip, with 
special importance for the last atom of tip. Also the parame- 
ter W2 could be modeled more accurately by introducing the 
spatial dependence between the tip and the graphene surface. 
Density functional calculations can be used to model the tip 
in realistic way. 

In our analysis we have consider that the electrons can only 
tunnel from the tip to the atom underneath, but when the tip 
is between two atoms the tunnel will take place to more than 
one atom. However, in this case the amplitude W2 is strongly 
reduced and the tunneling current will show a minimum. In 
fact it can be shown that the self energy has in this case the 
form 



AA 



(67) 



assuming the tip exactly in between the A and B atoms, and 
therefore the value of the current will be controlled mainly by 
the value of W2 ■ 

As we have seen in the calculation of the density of states, 
the case of weak coupling leads to the appearance of reso- 
nances at large energy values. The scan of the bias in our 
calculations did not reach these energies, but if it does a clear 
signature will be seen in the tunneling current. Also, for fixed 
bias, moving the tip away form the impurity will show oscil- 
lations in the STM current. 

One ingredient not included in our model is the effect of 
on-site energies at the impurities.^ As a consequence, a nat- 
ural question is whether our results are strongly modified if 
this effect is included. One should note that the values of the 
on-site energy and the hopping are correlated with each other 
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in the case of boron and nitrogen impurities. When the hop- 
ping in enhanced the on-site energy is positive relatively to 
that of the carbon atoms. In this case, the results of our calcu- 
lations in the present work show no qualitative changes from 
the more general model, as more detailed calculations show. 
In the case of reduced hopping, the same resonances we see in 
the local density of states of our work are also present in the 
more general approach, but their position in energy changes 
and an asymmetry of those resonances relatively to the Dirac 
point develops. One should also note that the different chem- 
ical nature of the atoms making the tip of the microscope also 
induces an asymmetry in the current, even when the LDOS 
does not show such asymmetry. Therefore, a careful study is 
needed to disentangle the effects from the tip and from a finite 
on-site energy at the impurity atoms. 

Another relevant question is whether the the resonances see 
in the density of state will broaden so much when the con- 
centration of impurities increases leaving no trace of them 
in the LDOS. Our calculations naturally refer to the diluted 
limit, where the distance between impurities is large. When 
the concentration increases the resonances will broaden, but 
their finger prints still remains in the LDOS."*** In Ref. 48, Fig. 
13(a) shows that the resonances remain well defined even for 
concentration of impurities up to 10%. It is not possible to 
give a characteristic length at the Dirac point such that above 
it the impurities could be considered to act isolated. This is 
so because Fermi momentum is zero. One was therefore to 
rely on numerical calculations with a varying concentration 
of impurities as in Ref. l48l 
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APPENDIX A: FREE PROPAGATORS AND USEFUL 
RELATIONS 



where Ac = 3%/3ao/2 is the area of the unit cell. With the 
above definition p{E) is finite in the energy range Q < E < 
3t. In addition we also have the useful relation 



The propagators for the perfect lattice are {h = 1) 



G%{uJn,k,p) : 



(Zc^„)2„i2|0(fc)|2 

4.p0*(fc) 

(*u;„)2-t2|0(fc)|2 

(zc^„)2-t2|0(fc)|2 



(Al) 

(A2) 

(A3) 

(A4) 
(A5) 



^0(fc)G%K,fe) = ^0*(fc)G%(a;„,fc) 



— ^AA\y>n) ■ 



(A9) 



Similar equations hold for G^^. 



* K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, 666 (2004). 

S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 



10 



^ K. S. Novoselov, D. Jiang, T. Booth, V.V. Khotkevich, S. M. Mo- 

rozov, A. K. Geim, PNAS 102, 10451 (2005). 
^ A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, 

and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009). 

N. M. R. Peres, F. Guinea, and A. H. Castro Neto, Phys. Rev. B 

73, 125411 (2006). 
' V. M. Pereira, F. Guinea, J. M. B. Lopes dos Santos, N. M. 

R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 96, (036801) 

(2006). 

^ Yu. V. Skrypnyk and V. M. Loktev, Phys. Rev. B 73, 241402 
(2006). 

V. V. Cheianov and V. I. Falko, Phys. Rev. Lett. 97, 226801 

(2006) . 

^ Cristina Bena, Phys. Rev. Lett. 100, 076601 (2008). 
' M. A. H. Vozmediano, M. P. Lopez-Sancho, T. Stauber, F. Guinea, 
Phys. Rev. B 72, 155121 (2005). 

A. Cortijo and M. A. H. Vozmediano, Nucl. Phys. B763, 293 

(2007) . 

" O. V. Yazyev, I. Tavernelli, U. Rothlisberger, and L. Helm, Phys. 
Rev. B 75, 115418 (2007). 

O. V. Yazyev and L. Helm, Phys. Rev. B 75, 125408 (2007). 
" N. M. R. Peres, A. H. Castro Neto, and F. Guinea, Phys. Rev. B 
73, 195411 (2006). 

E. R. Mucciolo, A. H. Castro Neto, and C. H. Lewenkopf, "Con- 
ductance quantization and transport gap in disordered graphene 
nanoribbons" , arXiv:0806.3777 1'l. 

E. Stolyarova, K. T. Rim, S. Ryu, J. Maultzsch, P Kim, L. E. Brus, 
T. F. Heinz, M. S. Hybertsen, and G. W. Flynn, Proc. Natl. Acad. 
Sci. USA 104, 9209 (2007). 

M. I. Katsnelson and A. K. Geim, Phil. Trans. R. Soc. A 366, 195 

(2008) . 

A. Fasolino, J. H. Los, and M. I. Katsnelson, "Intrinsic ripples in 
graphene". Nature Materials 6, 858 (2007). 

J. Martin, N. Akerman, G. Ulbricht, T. Lohmann, J. H. Smet, K. 
von Klitzing and A. Yacoby, Nature Phys. 4, 144 (2008). 

D. C. Elias, R. R. Nair, T. M. G. Mohiuddin, S. V. Morozov, P 
Blake, M. P Halsall, A. C. Ferrari, D. W. Boukhvalov, M. I. Kat- 
snelson, A. K. Geim, and K. S. Novoselov, "Control ofgraphene's 
properties by reversible hydrogenation" , arXiv:08 10.4706] 

E. J. Mele and J. J. Ritsko, Phys. Rev. B 24, 1000 (1981). 

^' M. Endo, T. Hayashi, S.-H. Hong, T. Enoki, and M. S. Dressel- 
haus, J. Appl. Phys. 90, 5670 (2001). 

O. Stephan, P M. Ajayan, C. Colliex, P. Redlich, J. M. Lambert, 
P Bernier, and P Lefin, Science 266, 1683 (1994). 
T. B. Martins, R. H. Miwa, A. J. R. da Silva, and A. Fazzio, Phys. 
Rev. Lett. 98, 196803 (2007). 

N. M. R. Peres, F. D. Klironomos, S.-W. Tsai, J. R. Santos, J. M. 

B. Lopes dos Santos, and A. H. Castro Neto, Europhys. Lett. 80, 
67007 (2007). 

G. M. Rutter, J. N. Grain, N. P Guisinger, T. Li, P N. First, and J. 
A. Stroscio, Science 317, 219 (2007). 

P Mallet, F. Varchon, C. Naud, L. Magaud, C. Berger, and J.-Y. 
Veuillen, Phys. Rev. B 76, 041403(R) (2007). 



V. W. Brar, Y. Zhang, Y. Yayon, T. Ohta, J. L. McChesney, A. 
Bostwick, E. Rotenberg, K. Horn, and M. F. Crommie, Appl. 
Phys. Lett. 91, 122102 (2007). 

M. Ishigami, J. H. Chen, W. G. Cullen, M. S. Fuhrer, and E. D. 
Williams, Nano Lett. 7, 1643 (2007). 

V. Geringer, M. Liebmann, T. Echtermeyer, S. Runte, M. Schmidt, 
R. Rckamp, M. Lemme, and M. Morgenstern, "Intrinsic and ex- 
trinsic corrugation of monolayer graphene deposited on Si02", 
arXiv:0806.1028 (2008). 
^° Y. Zhang, V. W. Brar, F. Wang, C. Girit, Y. Yayon, M. Panlasigui, 
A. Zetl, and M. F. Crommie, Nature Phys. 4, 627 (2008). 
A. Deshpande, W. Bao, F. Miao, C. N. Lau, and B. J. LeRoy, 
"Spatially resolved spectroscopy of monolayer graphene on 
Si02", arXiv:0812.1073 (2008). 

H. Amara, S. Latil, V. Meunier, Ph. Lambin, and J.-C. Charlier, 
Phys. Rev. B 76, 115423 (2007). 
" T. O. Wehling, A. V. Balatsky, M. I. Katsnelson, A. I. Lichtenstein, 
K. Schamberg, and R. Wiesendanger, Phys. Rev. B 75, 125425 
(2007). 

Yu. V. Skrypnyk and V. M. Loktev, Phys. Rev. B 75, 245401 

(2007) . 

V. Mujica, M. Kemp, and M. A. Ratner, J. Chem. Phys. 101, 6849 
(1994). 

Taturo Fukuda, Huseyin Oymak, and Jongbae Hong Phys. Rev. B 
75, 195428 (2007); idem, J. Phys. Condens. Matter 20, 055207 

(2008) . 

" T. Stauber, N. M. R. Peres, and A. K. Geim, Phys. Rev. B 78, 
085432 (2008). 

K. S. Dy, Shi-Yu Wu, and T. Spratlin, Phys. Rev. B 20, 4237 
(1979). 

John Tomfohr and Otto F Sankey, J. Chem. Phys. 120, 1542 
(2004). 

N. M. R. Peres, T. Stauber, and J. M. B. Lopes dos Santos, Phys. 
Rev. B 79, 035107 (2009). 

C. Caroli, R. Combescot, P. Nozieres, and D. Saint-James, J. Phys. 
C:Solid St. Phys., 4, 916 (1971); idem, ibidem, 4, 2598 (1971); 
idem, ibidem, 5, 21 (1972). 

David K. Ferry and Stephen M. Goodnick, Transport in Nanos- 
tructures, (Cambridge: Cambridge University Press, 2001). 
'^^ Hartmut Haug and Antti-Pekka Jauho, Quantum Kinetics in 
Transport and Optics of Semiconductors, 2'^'^, (Springer: Berlin, 
2008). 

Abhishek Dhar and Diptiman Sen, Phys. Rev. B 73, 085119 
(2006). 

Carlo Beenakker and Christian Schonenberger, Physics Today 5, 
37 (2003). 

Ya. M. Blanter and M. Buttiker, Physics Reports 336, 1 (2000). 
"* '' Sylvain Latil, Stephan Roche, Didier Mayou, and Jean-Christophe 
Charlier, Phys. Rev. Lett. 92, 256805 (2004). 
Vitor M. Pereira, J. M. B. Lopes dos Santos, and A. H. Castro 
Neto, Phys. Rev. B 77 115109 (2008). 



